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I.  INTRODUCTION 

Given  • nonlinear  parachute  gliding  system,  control  can  he  effected  by  a 
servomotor  pulling  on  the  shroud  lines  of  the  parachute  which  causes  a banked 
turn  of  the  parachute  and  a change  in  the  direction  of  flight.  Under  the  condition 
that  the  magnitude  of  the  wind  Telocity  rector  w be  less  than  that  of  the  para- 
chute horizontal  reloci ty  rector  v (relatire  to  air);  i.e.,  ||w||  < ||v||,  the 
parachute  possesses  a wind  penetrating  capability  and  the  potential  of  reaching 
the  target  under  arbitrary  wind  directions.  It  is  assumed  that  the  initial 
altitude,  and  thus  the  gliding  period,  is  properly  chosen  in  accordance  with  the 
wind  angle  0 and  the  ratio  ||w||/||r||  such  that  a solution  exists.  Also,  a 
uniform  wind,  constant  in  both  magnitude  and  direction,  is  assumed  throughout 
this  analysis. 

It  is  desired  to  compute  a control  law  which  minimizes  the  control  effort 
and  the  terminal  error  from  the  target.  The  usual  linearization  analysis  is 
insufficient  to  approximate  the  given  system  due  to  the  highly  transcendental 
nonlinearities.  Thus,  a numerical  solution  to  this  problem  is  developed  and 
intended  to  give  some  light  to  the  stochastic  wind  case  study.  The  numerical 
algorithm  derived  by  Martensson  [4]  in  solving  the  Hamilton- Jacobi-Bellman 
partial  differential  equation  with  the  differential  dynamic  programming  (DDP) 
principle  is  applied  to  solve  the  given  nonlinear  gliding  problem.  The 
theoretical  derivation  of  this  algorithm  is  given  in  the  Appendix. 

II.  PROBLEM  FORMULATION 

A.  System  Equations  and  Transformation 

Assuming  a constant  wind  w with  angle  6 in  the  horizontal  plane,  a con- 
stant rate  of  de«.r*r.t  v,  and  an  initial  altitude  hQ  at  launch  time,  the  equations 
of  motion  governing  the  parachute  can  be  expressed  in  the  horizontal  plane  as 
follows:  (Pearson  [5]) 
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P ■ * ♦ * 0 ( t t T = —•  (] 

p:  position  vector 

v:  horizontal  velocity  cf  parachute  relative  to  air. 

This  is  depicted  in  Fig.  1.  The  velocity  vector  v(t)  is  assumed  to  have  con- 
stant magnitude  a.  Thus  v can  be  represented  by 


v^t)  = a cosCw(t)] 
v2(t)  * a sin[(ii(t)] 


where  the  velocity  angle  u(t)  is  related  to  the  bank  angle  4 of  the  parachute 
via  i:he  well-known  relation 


u(t ) = £ tan4(t) 


g:  gravity  acceleration. 


Since  the  bank  angle  4 can  be  directly  manipulated  by  changes  in  the 
servomotor  connecting  the  shroud  lines , we  can  rewrite  eq.  (1)  - (3)  as 


d 

dt  P1 


dt  p2 


= a cosu  + w. 


= a sinm  + w„ 


in  which  u is  regarded  as  the  control  variable.  Let  p(t  ),  w(t  ),  and  w(t  ) 

o o o 

be  given  at  some  initial  time  t in  the  interval  0 < t < T.  A performance 

o o 

index  which  takes  into  account  several  desirable  features  of  this  problem  is 


P(u)  = 5-  |jp(T)|| 


|2  + 5-  P rtt(T)-e-*]2  + 2TFTJ  f u2(t)dt 

0 \ 
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where  q and  are  non-negative  weighting  parameters.  In  eq.  (5)  the  first 
ten,  1/2  ||p(T)j|2,  reflects  the  desirability  of  Minimizing  the  Euclidean  dis- 
tance from  the  target  at  the  predetermined  terminal  time  T.  The  second  term 
reflects  the  desirability  of  hawing  the  parachute  point  upwind  at  the  terminal 
time  in  order  to  reduce  the  net  horizontal  velocity  and  thus  the  impact  at  touch- 
down. The  last  term  reflects  the  cost  of  control  effort  in  terms  of  the  'average 

power'  spent  over  the  interval  t ( t f T. 

o 

A time-varying  transformation  of  the  origin  could  be  made  according  to 


y s p + (T-t)w  tQ  < t ( T (6) 

Then  minimizing  ||p(T)||  is  equivalent  to  minimizing  ||y(T)||.  Moreover,  the 
independent  variable  t could  be  transformed  via 


Define  a set  of  new  variables 


V x2*  x3  by 


(7) 


. yl  y2 

X1  a(T-t  ) * x2  ' a(T-t  ) ' x3  " “ 
o o 


(8) 


The  system  equations  after  transformation  become 

*1  s 008  x3 

*2  = sin  x3  0 ( t ( 1 (9) 

* 

s (T-t  )u  = u 
o o 

with  initial  conditions 

V0)  * + <T-to)"i] 

o 

V°>  * OT3TT  ^2^0^  + <T'W 

O 

x,(0)  = w(t  ) 

O O 
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Here  denotes  differentiation  with  respect  to  normalised  time  t.  The  per- 
formance index  in  terms  of  these  new  variables  can  be  rewritten  as 


P(u)  * x2(l)  + xj(l)  + Q1[x3(l)-e-ir]2  + Q2  f1  u?(T)dT 

* A 


(10) 


V-2 


a (T-t  Y 
o 


In  vector  form  we  have 


°2  5nr  4 

a iT-t  r 
o 


x = A(x)  ♦ Bu 


where 


P(u) 


[A°(x)  ♦ B°(u,t)]dt 


COS  Xg 

k a 
0 

A(x)  = | 

sin  x3 

, B = 

0 

0 

1 

A°(x)  = x2(l)  * x2(l)  ♦ Q^XgCD-e-x]2 
B°(u,t)  = 02u2(t)  . 


(11) 

(12) 


B.  Existence  of  Optimal  Control 

According  to  an  existence  theorem  by  Lee  and  Markus  [3],  we  can  show  the 
existence  of  an  optimal  control  to  our  problem. 

Theorem  2 Consider  the  following  process  in  Rn: 

* = A(x,t)  + B(x,t)u  0 ? t f T 
A performance  index 

P(u)  = [ [A°(x,t)  + B°(u,t)]dt 

is  to  be  minimized. 

Assume: 

(i)  A,  B,  A0 , B° , ||,  |i  are  continuous  for  all  x e Rn,  u e Rm,  t e R 
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(ii)  A°(x,t)  i 0 V(x,t)  c Rn+i 

(iii)  B°(u,t)  i «|u|^  for  some  constant  w > 0,  p > 1 

(iv)  For  each  fixed  t,  B°(u,t)  is  convex  in  u 

(v)  |x(t)|  ( 8(|u|1),  where  8(0  is  a monotonic  increasing  function. 

(|»|^  expresses  L^-norm). 

Let  the  set  of  admissible  controls  be  L^[0,T]  such  that  the  response  x(t ) 
initiating  at  xq  yields  a finite  P(u).  Then  there  exists  an  optimal  control  u* 
which  minimizes  P(u). 

Proof  Refer  to  [3]. 

In  view  of  eqs.  (11)  and  (12),  conditions  (i)  thru  (iv)  are  evidently 
satisfied.  Let 


|x(t)|  a l | x . ( t ) I 0(tJl 

i=l  1 


Since 


*,(t)  * x1(o)  + cos[xg(s)]ds 


it  follows  from  the  triangle  inequality  that 


«x(t) | * lx1(o)|  + | |cos[x3(s)]|ds  * ix1(o)|  + 1. 


Similarly, 


|x2(t)|  * |x2(o)|  + 1 
and  ^ 

|Xg(t)|  ( 1 Xg(o) | + [ |u(s)jdS  = | Xg ( O ) | + |u(t)|.  . 

J o 

Substituting  these  inequalities  into  eq.  (13), 

3 

| x(t ) ! $2  + 1 x,(o)  + |u(t)L  5 B(|u(t)|.)  . 

1 A 1 1 

Obviously,  8(0  is  monotonically  increasing  in  its  argument;  thus,  condition  (v) 
is  also  fulfilled.  By  Theorem  2.1  there  exists  an  optimal  control 
u*  e L.[o,l]  which  minimizes  the  performance  index  P(u). 
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C.  Necessary  Conditions  for  Terminal  Constraint  Problem 
A more  appealing  problem  is  to  incorporate  the  terminal  constraints  into 
the  original  optimization  problem,  i.e.,  to  consider  the  following  modified 
performance  index: 

P(u)  = x2(l)  + x2(l)  + Q^XgdM-ir]2  + b^d)  + b2x3(l) 


+ b3Cx,d)-e-»]  ♦ Q2  [ u2(s)ds 

* ft 


(14) 


b.'s  are  appropriate  Lagrange  multipliers. 

In  this  constrained  optimization  formulation,  we  require  that  the  parachute 
be  driven  to  the  target  and  directed  opposite  to  the  wind  at  touchdown,  while  the 
average  energy  spent  is  minimized.  Recall  the  integral  form  of  die  system 
equations 

x1<t)  = x1(o)  + j cos[x3(s)]ds 


x2(t)  = x2(o)  + f sin[x3(s)]ds 
; o 


0{t(l 


(15) 


By  the  Schwartz  inequality. 


Ix^t)  - x1(o)|2  ( (f  |cos[x3(s)]|ds}2  * f cos2[x3(s)]ds 

; o 'o 

)|2  * f 

• rs 


ix2(t)  - x2(o 


sin  [x3(s)]ds  . 


This  implies 


1 x^ ( 1 ) - X1(o)|2  + I x2 ( 1 ) - X2(o)j2  ( 1 . 

Hence,  if  a control  exists  which  drives  x1  and  x2  to  the  origin  at  the  terminal 
time,  it  is  necessary  that 

|x1(o)|2  + |x2(o)|2  f 1 

In  other  words,  the  initial  conditions  of  x^  and  x2  must  lie  within  the  unit 
circle  on  the  horizontal  normalized  plane. 
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One  other  observation  is  that  if  we  negate  the  trajectory  of  xg  by 
applying  a negative  control  and  a negative  initial  condition  on  x3,  the  system 
equations  become 


-Xg(t)  = -XgCC^  + 

f C-u(s)]ds 

^0 

-Xj(t)  = -x2(o)  + 

[ sin[-x,(s)!Jd8 

'o 

x^(t)  * x^o)  + 

| cos[-x3(s)]ds 

Due  to  the  symmetric  property  of  the  cosine  function,  we  conclude  that  it  is 
sufficient  to  consider  the  upper  half  unit  circle  as  the  wrrking  zone  of  our 
problem.  This  is  shown  in  Fig.  2.  The  negation  of  an  op  . :.«1  control  u*  for  a 
set  of  initial  conditions  (x^o),  x2<o),  Xgfo))  would  be  an  optimal  control  for 
the  set  of  initial  conditions  (x^o),  -x^c),  -Xg(o)). 


III.  COMPUTER  SIMULATION 

The  purpose  of  this  simulation  work  is  to  study  the  practicality  of  comput- 
ing an  optimal  control  on-line  with  the  DDP  algorithm.  Given  a uniform  wind  with 
known  magnitude  and  direction,  a constant  descent  rate,  initial  altitude  and 
direction  of  the  parachute  velocity  vector  relative  to  air,  selected  points  in 
the  upper  unit  circle  are  chosen  with  various  launching  angle  u(o)  to  initiate  a 
search  for  the  optimal  control  as  a function  of  time.  Since  the  initial  guess 
of  the  nominal  control  and  multipliers  plays  a crucial  role  in  the  convergence , 
in  order  to  make  a proper  correction  of  an  initial  guess  promptly  after  the  first 
one  fails , the  entire  simulations  have  been  executed  under  a t ime-sharing  system 
CP/CMS  of  IBM/360. 
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A.  Simulation  Procedure 


In  accordance  with  the  DOP  algorithm  deacribed  in  the  Appendix,  an  outline 
of  the  simulation  procedure  ia  as  follows: 

1.  Guess  a nominal  control  u(t),  0 ( t ( 1,  and  integrate  the  system  equations  to 

obtain  x(t).  Store  u(t),  x(t).  Guess  a set  of  nominal  multipliers  b,  and 

compute  the  corresponding  nominal  cost  V(x  ,b,0). 

o 

2.  Compute  boundary  conditions  for  "a",  V , and  V , then  integrate  4,  V , V 

X XX  X XX 

backwards  from  t * 1 to  t * 0,  while  minimizing  H(x,u,Vx,t).  Compute 
and  store  the  minimizing  control  u*(t),  $ , and  a(x,b,t). 

3.  Compare  |a(xo,b,0)|  with  a specified  small  quantity  n . If  |a|  < i^,  the 
predicted  change  in  the  performance  index  is  small;  (u,x)  thus  is  considered 
as  the  optimal  solution  of  a constraint-free  problem.  Otherwise,  go  to  4. 

In  addition,  if  | |ip(x(l) ,1) 1 1 < n2»  where  n2  is  a specified  allowable 
terminal  error,  then  (u,x,b)  is  considered  as  the  optimal  solution  of  the 
terminal  constraint  problem.  Stop. 

4.  Apply  the  modified  control  u = u*  + B^x-x)  from  t * 0 to  t = 1.  If  the 

reduction  in  cost,  V-V,  is  large  enough  compared  with  the  predicted  change  in 
V-V 

cost;  e.g.,  yjj-  > c (=  0.5),  c is  an  empirical  factor,  then  let  (u,x)  be  the 
new  nominal  solution  and  go  to  2.  Otherwise,  ’step-size  adjustment'  technique 
is  required.  (Jacobson  and  Mayne  [2]) 

5.  Modify  the  multipliers  by  6b  so  that  F + (5  i 6b)'ip  + f Lds  attains  a minimum 

■*0 

corresponding  to  ip  * 0.  Then  compute  boundary  conditions  for  V^,  V^,  V^. 
Integrate  V^,  V^,  Vbb  backward  from  t = 1 to  t = 0.  Compute  B2,  and  store  it. 

6.  To  compute  6b,  notice  that  for  the  optimal  solution 

V,(x  ,b*,0)  = <p'(x*(l),l)  = 0 where  b*  = b + 6b  . 

J5  O 

Expand  V^  to  first  order  about  b and  assume  *s  nonsin£ular» 

we  have 
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7. 


*■>  * -Vbk<1'o’®'0)Vb<xo*^'0' 

Apply  the  new  control  u * u*  + B^ix  ♦ $26b  to  the  system.  Store  u and  x, 
then  check  whether  6b  is  acceptable  or  not  according  to 
(1)  ||t(x(l),l)||  - ||*<x(l),l)||  > 0 

V(x  ,b*,0)  - V(x  .6,0) 

(11)  Y,  > > Y, 

V(x  ,b*,0)  - V(x  .6,0)  A 
o o 

where  V Is  the  new  nominal  cost,  y1  and  y2  are  suitably  chosen  (e.g., 

0 < Yj  < 1*  Y2  > D*  If  satisfied,  (u,x)  Is  a new  improved  nominal  solution 

with  cost  V.  Return  to  2.  If  not,  go  to  8. 

8.  Choose  (6b)  ^ * 4 6b  and  go  to  7.  If  no  correction  has  been  attained  after 
new  z 

several  reductions  of  6b,  then  (ii)  of  7 is  released.  The  only  demand  on 
6b  is  to  reduce  | U|  | . 

B.  Simulation  Results 

In  Section  II. B we  have  shown  that  it  is  sufficient  to  consider  the  upper 
half  unit  circle  as  a working  zone.  Twelve  rays  of  angles  0°,  14°,  30°,  45°, 
68.2°,  90°,  105°,  120°,  135°,  150°,  165°  and  180°  partition  the  upper-half  unit 
circle  into  eleven  sectors.  Five  points  are  chosen  along  each  ray  as  initial 
launching  points.  They  are  located  at  0.1,  0.35,  0.6,  0.8  and  0.95  of  the  unit 
radius  and  expressed  as  A,  B,  C,  D,  E,  respectively.  It  is  assumed  that  the 
complete  feasible  region  is  so  smooth  that  the  distribution  pattern  over  these 
sixty  points  is  enough  to  represent  the  upper  half  unit  circle.  Table  1 lists 
the  coordinates  of  these  points. 

Without  loss  of  generality,  an  easterly  wind  (along  the  negative  x^-axis) 
is  considered  of  magnitude  20  ft/sec.  The  magnitude  of  the  parachute  velocity 
relative  to  air  is  30  ft/sec.  An  error  criterion  of  0.015  is  considered  in  the 
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iterative  process,  which  allows  maximum  deviation  from  the  target  at  touchdown  of 
no  more  than  45  ft.  and  deviation  from  the  up-wind  direction  of  no  more  than  8.5°. 
Weighting  factors  Q1  and  Q2  are  1,0  and  0.1,  respectively.  The  gliding  time  is 
100  sec.;  however,  only  normalized  time  appears  in  the  simulation  process.  The 
unit  time  interval  is  once  divided  into  100  subintervals  when  an  integration 
routine  is  applied,  although  the  DDP  algorithr  is  originally  devised  to  equip  the 
integration  routine  with  500  sub intervals . By  taking  this  subdivision  it  is  found 
that  only  one-fourth  of  computer  storage  (128K  bytes)  and  one-third  of  cost  are 
required  in  addition  to  a loss  of  accuracy  within  an  order  of  1%  compared  to  that 
of  '500'  subdivision. 

Based  on  the  information  available  before  the  search,  different  initial 
guesses  of  nominal  control  and  multipliers  are  used.  Among  these,  constant 
control  with  null  multipliers  as  well  as  piecewise  constant  control  with  null 
multipliers  are  interchangeably  chosen  for  every  initial  setting  of  launching 
angle.  The  value  of  the  control  is  determined  such  that  the  launching  angle  is 
driven  opposite  to  the  wind  at  the  terminal  time  by  this  control.  In  addition, 
for  each  fixed  radial  point  consecutive  changes  of  launching  angle  are  made  to 
form  a complete  data  set.  At  each  step  change  of  these  angles,  a previous 
successful  optimal  control  and  multipliers  are  considered  as  current  initial 
guesses  of  the  nominal  control  and  multipliers.  In  most  instances,  as  one  would 
expect,  the  latter  combination  achieves  much  better  convergence  than  the  former 
two. 

It  is  observed  that  as  far  as  the  convergence  is  concerned,  there  is  a 
strong  dependence  on  the  way  launching  angle  is  measured.  More  specifically,  if 
we  define  a counterclockwise  measurement  to  be  positive,  then  a convergent  solu- 
tion is  achieved  more  frequently  for  a negative  launching  angle  than  for  its 
positive  counterpart  (i.e.,  + 2w).  This  can  be  seen  from  the  fact  that  the 

velocity  angle  is  directly  effected  by  control  and  the  wind  is  pointed  along  the 
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negative  x^-axia  which  requires  the  terminal  velocity  angle  to  be  zero.  A nega- 
tive measure  of  launching  angle  directs  the  trajectory  to  make  a counterclockwise 
turn  towards  the  wind  at  the  terminal  time.  For  instance,  consider  a test 
point  A in  the  first  quadrant  with  a launching  angle  -*  or  its  counterpart  it. 

We  can  plot  a tentative  optimal  trajectory  initiated  by  -t  as  the  dotted  line,  the 
trajectory  initiated  by  n as  the  circled  line  shewn  in  Fig.  3.  It  is  clear  from 
this  figure  that  the  dotted  trajectory  is  much  more  feasible  than  the  circled  one. 
In  this  case  a negative  launching  angle  with  a positive  control  which  gives  a 
counterclockwise  turn  towards  the  wind  has  a better  potential  to  reach  the 
origin  than  a trajectory  initiated  by  a positive  launching  angle.  Similar  situa- 
tions are  illustrated  at  test  points  B and  C.  This  geometrical  consideration 
determines  an  initial  guess  of  the  nominal  control  which  usually  results  in  a 
convergent  solution. 

Representative  optimal  controls  and  corresponding  trajectories  in  the 
normalized  plane  are  shown  in  Fig.  4 thru  Fig.  8.  Since  an  easterly  wind  is  con- 
sidered, the  velocity  angle  is  driven  to  zero  at  the  terminal  time.  This  is 
easily  seen  from  these  figures  by  noting  that  the  slope  of  x^-x^  trajectory  de- 
notes the  tangent  of  the  velocity  angle.  For  those  points  outside  80%  of  the  unit 
circle,  the  optimal  control  renders  most  of  its  effort  (in  terms  of  magnitude) 
to  bend  the  trajectory  at  the  beginning  towards  a natural  glide  (i.e.,  with  null 
control)  and  to  the  opposite  direction  of  wind  at  the  end  of  flight.  If  a bank 
angle  restriction  were  imposed  in  a practical  case,  which  means  the  control  is 
constrained  to  a certain  bounded  value,  then  a prolonged  gliding  period,  which  in 
turn  requires  a higher  launching  altitude,  has  to  be  used  in  place  of  the  100  sec. 
considered  here  to  scale  the  magnitude  of  this  control,  see  eq.  (9). 

It  is  also  interesting  to  notice  that  for  some  of  these  test  points,  the 
magnitude  of  launching  angle  is  even  raised  beyond  3u<j°  in  order  to  achieve  a 
convergent  solution.  This  trend  becomes  more  apparent  as  the  test  points  move 
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aiong  a radial  direct ion  towards  tha  origin.  Meanwhile,  the  optimal  control  pattern 
appears  to  exhibit  low  frequency  sinusoidal  characteristics.  In  other  words,  the 
control  effort  is  smoothed  throughout  the  whole  time  interval  rtther  i.han  sharply 
enforced  at  the  beginning  and  end  of  flight.  The  reason  is  that  the  closer  the 
test  point  to  the  origin,  the  more  effort  should  be  spent  to  counterbalance  the 
trade-off  between  the  fixed  gliding  time  and  distance  from  the  origin.  Thus  the 
optimal  trajectory  has  a somersault  tuning  effect;  i.e. , for  those  circumstances 
where  the  excess  flight  time  is  large  the  parachute  meanders  around  the  target. 

In  the  simulation  process,  it  is  felt  that  a good  choice  of  Lagrange 
multipliers  in  addition  to  the  nominal  control  usually  makes  a distinct  difference 
to  the  final  convergence.  The  self-adjustment  routines  provided  by  the  DDP 
algorithm  in  computing  the  improved  variation  of  multipliers  and  control  are 
based  on  a second  order  Taylor  expansion  about  the  optimal  solutions.  The  optimal 
multipliers  corresponding  to  the  optimal  sol'itions  depicted  in  Fig.  4 thru 
Fig.  8 are  given  in  Table  2.  From  this,  one  can  sketch  how  sensitive  these 
multipliers  are  to  the  initial  condition  changes. 

The  entire  simulation  work  is  summarized  in  Fig.  9.  Completely  and 
partially  feasible  regions  are  specified  individually  as  the  pear  shape  shaded 
area  and  funnel  shape  spots  outside  this  area.  This  map  provides  a corw»-  ^ent 
reference  far  the  pilot  to  drop  the  parachute  while  flying  into  the  completely 
feasible  region  under  uniform  wind  condition.  All  feasible  ranges  of  launching 
angle  are  tabulated  in  Table  3.  The  bank  angles  corresponding  to  the  magnitudes 
of  the  normalized  controls  are  shown  in  Fig.  10  for  various  (T-t0)  values.  A 
comparison  with  the  optimal  control  values  in  Figs.  4-8  indicates  that  bank 
angles  exceeding  30°  would  only  be  required  when  the  time  to  go  is  relatively 
short  — on  the  order  of  20  or  30  seconds. 

An  average  cost  of  $2.50  corresponding  to  a CPU  time  of  20  sec.  is  needed 
to  compute  a typical  solution. 
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IV. 


CONCLUSION  AND  REMARKS 


The  optimal  control  problem  for  a gliding  parachute  system  is  formulated 
in  two  ways.  One  is  to  minimize  the  performance  index  without  terminal  constraints, 
and  the  other  is  with  the  constraint.  We  have  shown  the  existence  of  an  optimal 
control  for  the  former  problem.  Nevertheless,  we  can  only  provide  some  convergent 
results  in  the  latter  case.  Two  reasons  may  explain  this.  First,  the  controlla- 
bility question  of  a nonlinear  system  of  the  type  considered  here  is  not  answered 
yet;  that  is,  the  existence  of  such  a control  that  drives  the  given  system  to 
the  origin  at  e prespecified  time  is  beyond  our  knowledge.  Second,  an  inherent 
deficiency  in  the  DDP  successive  approximation  scheme  arises  when  singular  matrix 
inversions  are  encountered. 

A completely  feasible  region  in  the  upper  half  unit  circle  is  computed  and 
specified.  In  th.'.s  region,  a convergent  solution  is  achieved  for  any  given  condi- 
tion. In  the  rest  of  the  upper  unit  circle,  for  different  initial  coordinates 
(x1,x2),  feasible  ranges  of  launching  angle  from  which  convergent  solutions  are 
obtained  have  also  been  specified. 

It  is  observed  that  the  optimal  control  function,  as  well  as  the  feasible 
range  of  launching  angle,  varies  considerably  along  both  the  radial  and  circular 
direction.  The  initial  guess  of  nominal  control  and  multipliers  thus  p’ays  an 
important  role  in  obtaining  convergence.  In  most  cases,  a geometrical  considera- 
tion helpful  to  determine  such  an  initial  guess. 

From  a practical  viewpoint,  further  investigation  is  needed  to  determine 
either  a more  efficient  algorithm  for  computing  optimal  controls  or  an  acceptable 
auboptimal  control  law  which  can  be  implemented  on-line. 
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APPENDIX 


Derivati.cn  of  a DDP  algorithm 

A second  order  differential  dynamic  programming  (DDP)  algorithm  ia  applied 
to  compute  the  optimal  control  for  selected  initial  conditions  in  the  upper  half 

unit  circle.  Her.  is  a brief  review  of  the  derivation  of  this  algorithm. 
(Martensson  [4]) 

Given  a dynamic  system 


* = f(x,  u,  t)  x(0)  = x . 

o 

and  the  performance  index 

( T 

u(U  J<U>  = u(?)  (r[x<T)*T:l  * { «x,u,t)dt} 

0<t<T  0<t<T  ' 


(16) 


(17) 


Subject  to  the  constraint  equation 


^Cx(T) ,T]  = 0 . (18 

One  way  to  manage  the  terminal  constraints  is  to  incorporate  * into  the  per- 
formance index  by  means  of  Lagrange  multipliers  (Bryson  and  Ho  tlj).  from  now  on 
W*  win  consider  ,e  modified  performance  index 

J(u)  = F[X(T),T]  + J L(x,u,t)dt  + b * Lx(T),T]  (19 

(')  means  taking  the  transpose  of  a vector  or  matrix.) 

De fine 


( T 

V(x,b,t)  = min  jF  + b if/  + f Lds 
u(s)  ( Jt 

t<s<T 


(20) 
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Assuming  that  V(x,b,t)  exists  and  is  twice  continuously  differentiable  with 
respect  to  x and  t for  all  tc[0,T],  then  V satisfies  the  Hamilton-Jacobi- 
Bellman  equation 


- It  1 mi"  {L  + vx  f} 

u v ' 


(21) 


Suppose  u,  x,  b is  a nominal  solution  neighboring  to  the  optimal  solution  u (t) 

_ mm  j'( 

= u(t)  + 5u(t),  x (t)  = x(t)  + 6x(t),  b = b + 6b  then  equation  (21)  becomes 


3V  ft  ft 
(x  ,b  ,t) 


min  (L(x  + 5x,u  + fiu,  t)  t V (x  + 5x,  b ♦ fib,  t)  f(x  + fix,  n + 6u,  t) 
fiu 

(22) 


it  it 

Now  assume  V(x  , b , t;  is  sufficiently  smooth  to  be  expanded  in  a second  order 

. ft  ft 

Taylor  expansion.  Then  we  can  approximate  V(x  , b , t)  with 

V(x“,  b",  t)  = V + Vx5x  + Vb5b  + <5x,  Vxfa5b>  + i <5x,  Vxx5x>  + j <5b,  Vbbfib> 


V (x'\  b*.  t)  = V + fix"  V ♦ fib"  V , 
x * x xx  xb 


(23) 


all  quantities  are  evaluated  at  x,  b,  t unless  otherwise  specified.  Let  a(x,  b,  t) 
be  the  difference  between  optimal  cost  at  x,  b,  t and  predicted  optimal  cost  at 
this  point,  i.e. 


a(x,  b,  t)  = V(x,  b , t)  - V(x,  b,  t) 


(24) 


Substituting  Eq.  (2C)-(24)  into  eq.  (22), 

3V._  3VU  3V 


■ {I?  + JT  6x  ♦ JT  5b  ♦ <«x,  fi’0>  + i <fix,  ^ 5x>  ♦ i <5b,  5b>! 


3V 

> 

3t 


3t 


min  {L  + < V + fix  V + fib  V , , f(x  + fix,  u + fiu,  t)>i 
6u  x xx  xd  j 


(25) 
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Because  V(x  , b , t)  is  approximated  by  a second  order  expansion,  the  following 
relations  between  total  and  partial  derivatives  holds 


~ (V  + a)  = •”■( V + a)  + f ( x , u,  t) 


d_ 

dt 


if*  1 “• l) 


d_ 

dt 


. 3V 

d „ _ xx 

dt  Vxx  " at 


3V, 

— * f “•  Vbx 


d_ 

dt 


(26) 


d_ 

dt 


All  V and  its  partial  derivatives  are  evaluated  at  x,  b,  t.  Eq.  (25)  and  (26) 
are  fundamental  equations  in  this  DDP  algorithm.  From  these  equations  we  can 
identify  all  those  partial  derivatives  of  V(x,  b,  t)  in  terms  of  the  Hamiltonian, 
which  we'll  introduce  next,  and  its  partial  derivatives  with  respect  to  u and  x. 
These  values  can  be  computed  simultaneously  when  we  integrate  Eq.  (26)  to  find 
V(x,  b,  t)  and  its  partial  derivatives  with  respect  to  x and  b. 

it 

Finally,  compute  V(x  , b , t)  according  to  Eq.  (23).  During  this  process, 

<5u  the  optimal  variation  which  should  be  added  to  the  minimizing  control  u 
(which  minimizes  Hamiltonian)  is  also  computed  in  terms  of  5x  and  6b. 

Define  the  Hamiltonian  to  be 


H(x,  u,  V , t)  = L(x,  u,  t)  + V f(x,  u,  t)  (27) 

A X 

By  this  definition,  if  we  let  6x  and  6b  be  zero  ir.  L'ne  right  hand  side  of  Eq.  (25), 
we  have  in  terms  of  H 
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(28) 


min  {H(x,  u + 6u,  V , t)} 

6u  X 

JL 

First  we  can  determine  the  optimal  variation  5u  that  minimizes  the  bracket 

. « - _ ft 

in  (28),  and  let  the  minimum  be  H(x,  u,  V , t)  where  u = u + 5u  . The 
necessary  conditions  for  this  are 


Hu(x,  u,  Vx  , t)  = 0 


and 


(29) 


Huu(x,  u,  Vx,  t)  > 0 (positive  definite) 

This  u would  be  the  optimal  solution  if  the  corresponding  trajectory  and 
multipliers  were  x and  b respectively.  However,  this  is  not  the  case  in  general. 
Therefore,  certain  correctiors on  u must  be  made  to  take  into  account  the 
difference  between  (x,  5)  and  the  optimal  one.  Let  these  differences  be  6x 
and  6b  respectively.  Reconsider  the  minimization  of  Hamiltonian  as 


min  (H(x  + 6x,  u + 5u,  V (x  + 5x,  b + 6b  , t),  t)} 
6u  x 


(30) 


Again,  necessary  conditions  for  a minimum  would  be 


Hu(x  + 5x,  u + 5u,  (x  t 5x,  5 + 5b,  t),  t)  = 0 


H^(x  + 5x,  u + 5u,  (x  + 5x,  b + 5b,  t),  t)>  0 


(31) 


In  order  to  determine  5u  in  terms  of  5x  and  5b,  we  expand  Eq.  (31)  to  first 
order  about  x,  u,  b,  t.  Then 


H + H 5u  + H 5x  + f (V  5x  + V . 5b)  = 0 
u uu  ax  u xx  xb 


From  Eq.  (29)  Hu  = 0,  hence 


H 5u  + (H  + f V ) 5x  + f V , 5b  = 0 
uu  ux  u xx  u xb 
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Assume  H (x,  u,  V , t)  is  nonsingular,  then 

UU  X 


6u  = g^  6x  + ^ 6b 


(32) 


where  g = - H _1  (H  + f V ) 
1 uu  ux  u xx 


B2  5 - "uu'1  fu"  V*b 


(33) 


Insert  6u  into  Eq.  (30)  and  expand  it  to  second  order  about  x,  u,  b,  t;  using 


the  fact  iMx,  u,  V , t)  = 0,  we  obtain 


* * 


H + (H  + f V ) 6x  + f V . 6b  + <6x,{(H  + f V ) g_  + g,  H g„ 

x xx  xb  * ux  u xx  2 1 uu  2 


* (fK  * W vxb>  sb>  * <Sb-  <2  B2  "uuB2  * vxb  f B2>  4b> 


♦ <Sx,  (i  Hxx  ♦ * i ♦ S*  (Hux  . fu  ”Vxx>)  ««> 


(34) 


whsre  all  quantities  are  evaluated  at  x,  u,  E,  t. 

The  series  expansion  (34)  is  now  ready  for  identification  with  the  left  hand 
side  of  Eq.  (25).  We  can  identify  those  partial  derivatives  by  equating  the 
coefficients  of  power  terms  in  6x  and  6b.  Finally,  after  combining  with  Eq.  (26) 
we  obtain  a set  of  ordinary  differential  equations  from  which  V,  V^,  Vfa,  V^, 

Vxx*  Vbb  31,6  inteF>rated» 

- = H - H(x,  u,  Vx  t) 


dV 

- = »x  + [f  - f(x,  u,  t)]  V 


XX 


dvb  _ . 

d r = Cf  • f(x»  u*  t)]  vxb 
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w 


dV 


xb 


dt 


= Cfx+W  V; 


xb 


dV. 


bb 

-nr-  = - 30  H e„ 
dt  2 uu  2 


dV 

xx 

dt 


" Hxx  + fx  Vxx  + vxx  fx  + Sl  Huu  S1  + ei  (HUX  + fu  Vxx) 


+ (H  + f v ) a 

ux  U XX;  B1 


with  boundary  conditions 


a(x(T),  b,  T)  = 0 

Vx(x(T),  b,  T)  = Fx(x(T),  T)  * b'  ^ (x(T) , T) 

Vb(x(T),  b,  T)  = i (x(T) , T) 

Vxb(x(T)»  T)  * *x  (x(T),  T) 

Vbb(x(T)’  T)  = 0 

Vxx(x(T)*  T>  = Fxx(^T>.  T)  * b * (x(T) , T) 

XX 

All  terms  on  the  right  hand  side  of  Eq.  (35)  are  evaluated  at  x,  G,  b, 
unless  other  wise  specified. 


” yu PS 


(35) 


(36) 


t 
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Table  1.  Selected  Raya  and  Radial  Points 


Table  3.  Feasible  Ranges  of  Launching  Angle 


Feasible  ranges  of  launching  angle  according  to  x.-axis  symmetry.  See  Section  II. C 
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FIG.  I PARACHUTE  DYNAMICS  IN  HORIZONTAL  PLANE. 


FIG.  2 UPPER  HALF  UNIT  CIRCLE. 
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NEGATIVE  MEASURE  * ^ WIND 


FIG.  3 TENTATIVE  OPTIMAL  TRAJECTORIES  FOR  DIFFERENT  MEASUREMENT 
OF  THE  LAUNCHING  ANGLE. 


4(a)  OPTIMAL  CONTROLS  VS  NORMALIZED 
TIME. 


x,  (0)  ■ 0.95 
x2(0)  * 0 


.2 


FIG.  4(b)  OPTIMAL  TRAJECTORIES  IN 
NORMALIZED  COORDINATES. 
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KIG.  5(b)  OPTIMAL  TRAJECTORIES  IN  NORMALIZED 
COORDINATES. 
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(-270*) 

(-600*) 


FIG.  6(0)  OPTIMAL  CONTROLS  VS.  NORMALIZED 
TIME. 


x, ( 0)  « 0 
x2(0)  « 0.1 


FIG. 6(b)  OPTIMAL  TRAJECTORIES  IN 
NORMALIZED  COORDINATES. 


COMPLETE  FEASIBLE 
REGION 


34 


FIG.  9 A SUMMARY  CHART  OF  FEASIBLE  REGION  IN  THE  UPPER  HALF  UNIT  CIRCLE 


<P  ( BANK  ANGLE) 


FIG.  10  BANK  ANGLE  v s NORMALIZED  CONTROL 


